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ABSTRACT 


Imaging  physical  model  data  provides  a  good  test  for  an  inversion  algorithm.  The 
physical  model  data  are  real  wave  fields  and  do  not  include  the  simplifications  of 
synthetic  data.  Also,  the  parameters  of  the  model  are  known  beforehand  so  that  it  is 
easy  to  determine  how  well  the  inversion  works.  Here,  inversion  is  a  true  amplitude 
Kirchhoff  depth  migration  in  the  sense  that  the  amplitude  of  the  imaged  reflections  is 
proportional  to  the  reflection  coefficient.  Each  shot  record  in  a  physical  model  data  set 
is  inverted  separately  with  a  common  shot,  prestack  inversion  routine  with  a  laterally 
and  depth  variable  velocity  function.  Each  shot  record  inversion  forms  a  partial  image 
of  the  subsurface.  The  results  are  then  stacked  to  form  a  full  image  of  the  subsurface. 
The  physical  model  data  set  is  inverted  twice.  For  the  second  inversion,  the  output 
trace  spacing  is  half  the  spacing  for  the  first  inversion  and  the  output  aperture  is  three 
times  wider  than  in  the  first  inversion.  In  both  cases,  the  background  velocity  field  is 
nearly  identical  to  the  actual  model.  This  tests  the  inversion  procedure  independent  of 
-ity  analysis.  Both  inversions  accurately  position  reflectors  in  the  model  but  each 
performs  better  on  different  portions  of  the  data.  With  a  larger  inversion  output  zone, 
steeper  events  are  imaged  better  but  the  increased  migration  “smile”  noise  obliterates 
some  deeper  events.  Both  inversions  are  superior  to  a  migration  of  the  data  performed 
by  Marathon  Oil  Company.  Velocities  and  densities  of  the  model  were  estimated  from 
the  inversion  amplitudes.  The  estimated  velocities  from  the  two  uppermost  reflectors 
agreed  within  ten  precent  of  the  true  velocities  but  the  estimates  degraded  with  the 
deeper  reflectors.  The  accuracy  of  the  density  estimates  could  not  be  determined 
because  the  true  densities  were  unavailable. 

The  physical  model  data  are  also  analyzed  for  doubly  mode-converted  reflections. 
The  data  are  dip  filtered  to  try  to  separate  any  mode  converted  reflections  from 
primary  compressional  wave  reflections.  Stacking  velocity  analysis  is  used  to  determine 
whether  the  dip-filtered  data  contain  any  mode  converted  reflections.  No  mode 
converted  energy  is  evident  in  the  data. 
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1.1  Cross  section  of  data  model. 

1.2  Schematic  of  recording  array  for  each  shot  in  data. 

1.3  Sample  sho»,  record  from  the  Marathon  data.  AGC  with  a  200- 
sample  (0.8  s)  window  has  been  applied. 

1.4  Ray  paths  followed  by  PSSP  mode  converted  waves. 

2.1  The  interfaces  of  the  constant  velocity  layers  of  the  background 
velocity  function  for  the  first  inversion. 

2.2  Inversion  of  the  shot  record  in  Figure  1.3.  Inversion  is  applied  to  the 
ungained  record.  The  inversion  shows  partial  images  of  the  reflectors. 

2.3  Stack  of  inversions  of  all  shot  records.  The  reflectors  are  completely 
imaged,  but  the  sawtooth  reflector  is  hard  to  see. 

2.4  Stacked  section  of  Figure  2.3  after  applying  AGC  with  a  200-sample 
(8000  feet)  window. 

2.5  The  envelope  of  all  reflectors  having  the  same  reflection  time.  It  is 
an  ellipse  with  the  source  and  receiver  at  the  foci.  Steep  reflected 
would  be  farther  than  flat  reflectors  from  the  source- receiver 
midpoint. 

2.6  Ray  tracing  from  a  subsurface  point  (x  =  13000  feet,  z  —  ',800  feet)  to 
equally  spaced  points  on  the  surface  using  the  first  background 
velocity  function. 
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2.7  The  smoothed  interfaces  of  the  constant  velocity  layers  of  the 
background  velocity  function  for  the  second  inversion. 

2.8  Ray  tracing  from  a  subsurface  point  (x  =  13000  feet,  z  =  7800  feet)  to 
equally  spaced  points  on  the  surface  using  the  second  background 
velocity  function. 

2.9  Second  inversion  of  the  shot  record  shown  in  Figure  1.3.  The 
impulse  response  smiles  cover  a  greater  portion  of  the  inversion  than 
the  images. 

2.10  Stack  of  second  shot  inversions  after  AGC.  Steeper  events  have  been 
imaged  better  than  in  Figure  2.4  but  there  is  more  smile  noise. 

2.11  Finite-difference,  time  migration  by  the  University  of  Houston.  The 
fourth  and  fifth  sawteeth  are  not  imaged. 

2.12  Prestack,  finite-difference,  depth  migration  by  Marathon  Oil 
Company.  The  fifth  sawtooth  is  absent. 

2.13  Inversion  (0)  of  the  shot  record  in  Figure  1.3  for  parameter 
estimation. 

2.14  Angle  of  incidence  of  rays  at  the  first  interface.  The  solid  line  is  the 
inverse  cosine  of  the  ratio  of  0  to  0} .  The  dashed  line  is 
tan-1  (x/5400),  the  theoretical  angle  if  the  first  reflector  is  a  flat 
plane  a  z  =  2700  feet. 


2.15  Plot  for  estimating  velocity  and  density  from  the  inversion 

amplitudes  of  the  first  reflector.  A  =  — sin  “Ojc“  and, 

Y  =  (1  —  i?)2/(l  -I-  R)  2  x  cos20/p  2 c2.  The  slope  of  the  best  fit  line  is 

the  reciprocal  of  the  square  of  the  density.  The  A-intercept  is  the 
reciprocal  of  the  square  of  the  velocity. 

2.16  Plot  for  estimating  velocity  and  density  from  the  inversion 

amplitudes  with  the  area  under  the  filter  selected  so  that  the  velocity 
is  15750  feet/s. 
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2.17  Plot  for  estimating  velocity  and  density  from  the  inversion 

amplitudes  of  the  second  reflector. 

2.18  Plot  for  estimating  velocity  and  density  from  the  inversion 

amplitudes  of  the  third  reflector. 

3.1  Slant  stack  of  a  physical  model  shot  record  from  Tatham  et  al. 
(1983). 

3.2  Slant  stack  of  the  (ungained)  shot  record  shown  in  Figure  1.3. 

3.3  Slant  stack  in  Figure  3.2  after  AGC  with  a  200-sample  (0.8  s) 
window. 

3.4  Inverse  slant  stack  of  Figure  3.2  for  p  >  11.8  x  10-6  s/ft.  This  is  the 
shot  record  in  Figure  1.3  after  dip  Filtering. 

3.5  Stacking  velocity  analysis  of  Figure  3.4.  A  PSSP  event  should  be  a 
peak  on  or  near  the  curved  line. 

3.6  Stacking  velocity  analysis  of  dip  filtered  CMP  gathers.  A  PSSP 
event  should  be  a  peak  on  or  near  the  curved  line. 

3.7  Constant  velocity  stack  with  v,  =  8000  feet/s.  AGC  has  been 
applied. 

3.8  P  to  S  transmission  coefficient  with  c  =  11750  ft/s,  cx  =  15750  ft/s, 
and  6!=0.3ci.  There  is  no  density  contrast.  More  conversion 
occurs  for  incidence  angles  greater  than  the  critical  angle  of  48° . 


3.9  S  to  P  transmission  coefficient  with  the  same  parameters  as  in  Figure 
3.8.  Most  conversion  occurs  for  emergence  angles  greater  than  the 
critical  angle  of  48° . 
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Symbol  Definition 

A  WKBJ  Green’s  function  amplitude 

0,  Inversion  operators 

b  Spatial  kernel  of  0 

c  Compressional  wave  velocity 

F  Amplitude  response  of  a  bandpass  filter 

7  Singular  function  of  a  surface 

p  Ray  parameter 

R  Angularly  dependent  reflection  coefficient 

a  Running  parameter  along  rays 

0  Incidence  angle 

t  Intercept  time 

T  Traveltime  function 

Us  Scattered  data 

x  Distance 


z 


Depth 


1.  INTRODUCTION 


A  physical-model  seismic  experiment  is  a  scaled  down  laboratory  simulation  of  a 
real  seismic  survey  and  thus  physical  model  data  are  an  analog  of  a  seismic  survey  over 
a  simplified,  known  earth.  The  data  are  collected  at  ultrasonic  frequencies  with 
distances  and  time  in  the  laboratory  experiment  scaled  so  that  the  distances,  time, 
frequencies,  and  medium  velocities  are  reasonable  for  real  seismic  data. 

These  data  are  useful  for  testing  and  comparing  seismic  data  imaging  techniques 
(migration/inversion).  Seismic  data  modeling  and  imaging  methods  are  based  on 
theory  that  incorporates  simplifying  assumptions  about  the  wave  field.  If  an  imaging 
procedure  is  based  on  the  same  theory  as  the  modeling  procedure,  the  imaging 
procedure  is  merely  the  inverse  of  the  modeling  procedure,  and  thus  while  the  imaging 
may  work  perfectly  on  synthetic  data  from  the  modeling,  it  might  not  work  well  on 
field  data.  Physical  model  data  do  not  introduce  this  problem;  the  data  are  real  wave 
fields  as  are  field  data.  The  model  data  contain  all  wave  effects,  including  head  waves, 
near-field  effects,  mode  conversions,  and  diffractions,  that  may  have  been  omitted  in 
the  synthetic  modeling.  Testing  imaging  techniques  on  model  data  gives  a  better 
indication  of  how  the  techniques  will  perform  on  field  data.  However,  because  the 
models  are  simpler  than  the  real  earth  and  the  physical  parameters  are  known,  it  is 
easy  to  verify  how  well  the  imaging  techniques  accomplish  what  their  intended  goal  of 
accurately  imaging  the  subsurface. 

Marathon  Oil  Company,  which  prompted  this  research,  has  generated  physical 
model  data  which  it  has  sent  to  contractors  to  evaluate  their  imaging  (migration) 
techniques.  The  true  velocities  are  known  beforehand,  so  they  can  be  used  as  the 
migration  background  velocity  function.  With  this  information,  the  migration 
techniques  can  be  compared  independent  of  the  errors  associated  with  imperfect 
knowledge  of  velocity.  Marathon  donated  a  physical  model  data  set  to  the  Center  for 
Wave  Phenomena  so  that  we  might  try  our  inversion  on  the  data.  The  model  is 
structurally  complicated  enough  to  warrant  prestack  inversion.  The  first  part  of  this 
thesis  shows  the  application  of  Dong’s  (1989)  common  shot,  c(x,z),  prestack  inversion 
routine  to  the  physical  model  data.  The  term  c(x,z)  means  that  the  input  velocity 
varies  laterally  and  vertically.  I  will  show  that  this  inversion  method  provides  an 
accurate  reflector  map  of  the  model  and  compares  favorably  to  migrations  of  the  same 
data  by  other  methods. 


1.1  Description  of  the  Physical  Model  Data 

The  data,  collected  at  the  Seismic  Acoustics  Laboratory  (SAL)  at  the  University 
of  Houston  for  Marathon  Oil  Company,  will  be  referred  to  as  the  Marathon  data.  The 
data  were  collected  over  a  block  model  in  a  water  tank  (Figure  1.1).  Geologically,  the 
model  represents  a  salt  ridge  over  a  rifted  basement.  The  top  layer  shown  in  Figure  1.1 
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was  water.  (Note  that  the  velocity  for  the  first  layer,  vp  —  11750  ft/s,  is  the  scaled 
velocity,  hence  the  difference  from  the  true  water  velocity.)  The  lower  layers,  the 
actual  block,  were  various  epoxy  resins.  The  high  velocity,  vp  =  22410  ft/s,  third  layer 
is  the  modeled  salt  ridge.  SAL  provided  the  scaled  compressional  wave  velocities,  vp.  I 
did  not  know  the  densities  or  shear  wave  velocities,  except  for  the  top  layer,  because  it 
is  water.  Water’s  density  is  1  g/cm3  and  its  shear  wave  velocity  is  zero.  The  model 
varies  only  slightly  in  the  out-of-plane  (i.e.  y)  direction  so  that  the  2.5-D  assumption 
used  later  for  deriving  the  inversion  operator  should  be  adequate  for  these  data.  The 
dimensions  in  Figure  1.1  are  labeled  with  feet  instead  of  scale  feet.  For  brevity,  I  will 
use  feet  and  seconds  instead  of  scale  feet  and  scale  seconds  throughout  this  thesis. 

Figure  1.2  depicts  the  recording  geometry  of  each  shot  record  from  the  Marathon 
data.  The  data  consist  of  291  shot  records.  For  each  shot  there  were  48  receivers  in  an 
end-on  spread  to  the  right  of  the  shot.  The  near  receiver  offset  was  800  ft  and  the 
receiver  spacing  was  80  ft,  so  the  far  receiver  offset  was  4560  ft.  The  shot  spacing  was 
also  80  ft,  with  the  first  shot  located  at  x  =  0  ft  and  the  last  was  at  x  —  23200  ft. 
Although  the  shots  and  receivers  were  at  the  depth  z  =  0  ft  shown  in  Figure  1.1,  this  is 
not  the  water  surface.  They  were  submerged  sufficiently  so  that  no  reflections  from  the 
water  surface  were  recorded.  For  each  shot,  two  seconds  of  data  were  recorded 
sampled  at  4  ms.  The  recording  filter  passed  frequencies  between  2  and  60  Hz. 

Figure  1.3  is  a  sample  shot  record.  Automatic  gain  control  (AGC)  has  been 
applied  for  the  display  to  aid  in  viewing  events.  The  shot  location  is  x  =  2000  ft,  and 
the  receiver  spread  extends  from  z  =  2800  to  x  =  6560  ft.  After  the  earliest  event,  the 
direct  wave,  the  first  curved  event  is  the  water  bottom  reflection.  Reflections  from  the 
second  and  third  interfaces  are  the  next  two  events,  and  the  strong  event  at  about  1.45 
s  is  a  reflection  from  the  model  bottom.  The  reflection  from  the  sawtooth  interface 
does  not  produce  an  easily  identifiable  event  on  this  record. 


1.2  2.5-D  Inversion  Theory 

In  this  section  I  will  briefly  outline  the  theory  used  for  inverting  the  Marathon 
data.  More  detailed  descriptions  of  the  theory  may  be  found  in  Dong  (1989),  Bleistein 
(1987),  Bleistein  et  al.  (1987),  Docherty  (1987),  Sullivan  and  Cohen  (1987),  and  Cohen 
et  al.  (1986).  Inversion  and  migration  have  the  same  goal:  imaging  reflectors  in  their 
proper  location.  The  philosophies  of  the  techniques,  however,  differ. 

Zero-offset  migration  takes  advantage  of  the  exploding  reflectors  model.  This 
model  has  two  aspects.  First,  the  earth  velocities  are  replaced  by  velocities  one  half  the 
true  value.  Second,  at  time  zero,  the  reflectors  “explode”  to  generate  an  upgoing  wave 
field,  which  is  recorded  at  the  surface.  These  exploding  reflector  data  are  an 
approximation  of  the  data  that  would  be  recorded  in  a  true  zero  offset  experiment.  To 
migrate  the  zero-offset  data,  the  data  are  extrapolated  downward  to  all  output  depth 
levels  using  the  half  velocity  of  the  exploding  reflectors  model.  This  is  a  transformation 
of  the  data  to  what  would  have  been  recorded  if  the  receivers  were  at  depth.  At  each 
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depth  level,  the  data  are  imaged  at  time  zero  (when  the  reflectors  explode).  This 
produces  a  migrated  image  of  the  subsurface.  For  prestack  migration,  the  receiver  and 
source  data  are  both  downward  extrapolated  and  the  data  axe  imaged  at  time  zero.  At 
time  zero,  a  receiver  coincident  with  a  source  will  show  a  reflection  if  the  receiver  and 
source  axe  on  a  reflector.  Different  types  of  migration,  such  as  Kirchhoff,  finite 
difference,  and  Stolt,  all  follow  this  downward  extrapolation  and  imaging  procedure. 
The  methods  differ  in  their  approximations  to  the  wave  Field  extrapolator  and  in  the 
domain  [t  —  xox  u—k)  in  which  extrapolation  is  applied.  More  details  on  migration  may 
be  found  in  Claerbout  (1985)  or  Yilmaz  (1987). 

The  inversion  operator  is  derived  from  the  solution  of  the  forward  problem.  An 
integral  representation  of  the  scattered  data  is  obtained  by  solving  the  acoustic  wave  or 
Helmholtz  equation.  Either  the  Kirchhoff  or  Born  approximation  is  used  in  solving  the 
forward  problem.  I  will  describe  the  method  of  deriving  the  inversion  operator  when 
the  Kirchhoff  approximation  is  used.  The  Green’s  functions  in  the  representation  axe 
replaced  by  the  WKBJ  approximation  because  seismic  data  are  high  frequency.  That 
is,  the  length  scales  of  interest  in  the  medium,  such  as  depth  to  the  reflector,  radius  of 
curvature  of  the  reflector,  etc.,  are  much  larger  than  the  nominal  wavelength  in  the 
data.  The  forward  scattered  data  then  have  the  representation 

*/»(w»£)  =  *w  f  dx'A(x',()  etwT{z  ,  (1-1) 


when  the  Kirchhoff  approximation  is  used.  In  equation  (1.1)  the  x'  integral  is  over  the 
reflecting  surface.  £  is  a  parameterization  of  the  source  and  receiver  coordinates  x,  and 
xr .  A  (x,  £)  is  a  specific  amplitude  term  defined  by  the  theory.  It  includes  spreading 
and  transmission  losses.  T(x,  £)  is  the  sum  of  traveltimes  from  the  subsurface  point  to 
the  source  and  the  subsurface  point  to  the  receiver: 

T(x,£)  =  t{x,x„)  +  t(x,xr)  . 


When  the  Born  approximation  is  used,  the  iu>  in  equation  (1.1)  is  replaced  by  w2,  the 
x'  integral  is  over  volume,  and  the  kernel  A  is  different. 

Equation  (1.1)  has  mathematical  form  similar  to  a  Fourier  transform.  From 
experience  with  Fourier  transforms  and  the  migration  concept  of  propagating  the 
source  and  receivers  backward,  an  educated  guess  would  be  that  the  phase  of  the 
inversion  operator  would  be  the  opposite  that  of  the  forward  solution.  The  three 
dimensional  inversion  formula  is 

0{x)  =  fd(b3D(x,t)fdvF(u)e-iuT<*-Uu,(v,Q  .  (1-2) 


In  equation  (1.2),  F(u)  is  a  bandpass  filter  to  recognize  the  bandlimited  nature  of 
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seismic  data.  F(u>)  is  the  product  of  the  amplitude  spectra  of  all  filters  applied  to  the 
data.  This  includes  the  source  wavelet,  receiver  effects,  recording  instrument  filters, 
and  any  processing  filters.  The  amplitude  factor,  63^(i,  £),  is  determined  by 
substituting  for  U9  from  equation  (1.1)  and  requiring  that  the  inversion  output,  /?(x), 
be  asymptotically  equal  to  the  reflection  coefficient  for  the  output  point,  x,  on  the 
reflector.  This  amplitude  factor  corrects  the  amplitude  of  the  data  for  divergence  and 
transmission  losses  and  is  different  for  different  sortings  of  the  data.  That  is,  for 
data  sorted  into  common-offset  gathers  is  not  the  same  as  63£>  for  common-shot 
gathers. 

The  inversion  operator,  defined  by  equation  (1.2),  is  essentially  a  Kirchhoff 
migration  operator  with  a  special  amplitude  weight.  Careful  attention  is  paid  to  fhe 
amplitude  in  deriving  the  forward  solution  and  the  inversion  operator  so  that  the 
amplitude  in  the  results  can  be  faithfully  diagnos  „ic  of  characteristics  of  the  medium. 

Equation  (1.2)  is  the  3-D  inversion  formula.  Most  seismic  data,  however,  are 
collected  in  2-D  lines.  For  a  single  seismic  line,  we  cannot  expect  to  be  able  to 
correctly  image  any  reflection  from  outside  the  plane  of  the  survey.  We,  therefore, 
assume  that  the  earth  velocity  varies  only  in  the  plane  directly  beneath  the  data  line 
and  has  no  out-of-plane  variation.  This  is  equivalent  to  assuming  that  the  seismic  line 
was  shot  along  the  dip  direction.  Seismic  sources,  however,  are  point  sources  and  this 
necessitates  the  use  of  3-D  Green’s  functions  that  take  three-dimensional  spreading  into 
account.  This  situation  with  2-D  data  and  3-D  (i.e.,  point)  sources  is  called 
2.5— dimensional. 

To  convert  the  3-D  inversion  formula,  equation  (1.2),  to  2.5-D,  the  integral  for  the 
out-of-plane  £— direction  is  evaluated  by  stationary  phase.  Stationary  phase  is  a 
method  for  asymptotic  evaluation  of  Fourier-like  integrals  Bleistein  (1984),  valid  under 
the  assumption  the  data  are  high  frequency.  After  the  stationary  phase  approximation 
has  been  applied,  we  obtain  the  2.5-D  inversion  formula: 

P(x)  -  J  d£b(x,£)f  du  F{u)\/iu  e~,wT(z'^  ,  (1-3) 

where  V 7u>  =  \/  j  u  |  e  ,s>%nM7r/4  .  This  factor  arises  in  the  stationary  phase  evaluation 
and  provides  the  proper  phase  for  the  inverted  data. 

As  stated  earlier,  the  inversion  asymptotically  gives  the  geometric  optics  reflection 
coefficient  when  the  output  point  is  on  the  reflector.  That  is, 

(5{x)  ~  R\x,d(x)}cos6(x)iD{x)  ,  (1.4) 


where  R\x,9(x)\  is  the  angularly  dependent  reflection  coefficient  and  7g(x)  is  the 
bandlimited  singular  function  of  the  reflector.  It  is  important  to  note  that  the  inversion 
output  gives  the  reflection  coefficient  at  only  one  angle  for  a  single  output  point.  This 
is  the  angle  for  which  a  specular  reflection  from  this  point  could  be  recorded  given  the 
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shot  and  receiver  geometry.  Therefore,  6  is  a  function  of  x  but  I  will  usually  not 
indicate  this  dependence.  That  is,  I  will  use  8  to  mean  6(x). 


The  singular  function  of  the  reflector  is  a  Dirac  delta  function  that  has  its  support 
on  a  surface  (Bleistein,  1984).  For  a  surface  S,  let  o  be  the  signed  normal  distance  from 
5.  Then,  the  singular  function  is  defined  as  7(1)  =  <5(cr).  On  the  surface, 


1b{x) 


b{ <*) 


tr  =  0 


=  JV(cj)  du  .  (1.5) 


Therefore,  the  inversion  output  for  a  point  on  the  reflector  is 

f3(z)  ~  R{x,8)  cos  8  J F(u)duj  ,  (1.6) 


and  for  a  point  not  on  a  reflector  it  decays  as  the  reciprocal  of  the  distance  to  the 
reflector. 

I  used  the  operator  defined  by  equation  (1.3)  to  invert  the  Marathon  data.  The 
implementation  of  the  operator  is  described  in  Chapter  2.  Further  details  may  be 
found  in  Dong  (1989). 


1.3  Why  Mode  Conversion? 

Shear  wave  (S-wave)  applications  to  exploration  have  increased  in  the  last  decade. 
Shear  wave  surveys  can  now  help  determine  anisotropy  or  can  be  used  for  direct 
hydrocarbon  detection  (Tatham  and  Stoffa,  1976;  Meissner  and  Hegazy,  1981; 
Robertson  and  Pritchett,  1985;  Neidell,  1985).  Marine  shear  wave  surveys  may  also 
have  applica  ons  in  seabed  mapping. 

On  land,  shear  waves  can  be  generated  and  recorded  directly  and  easily.  Since 
water  cannot  support  shear  waves,  marine  S-wave  surveys,  however,  are  difficult  to 
conduct.  Marine  shear  wave  data  can  be  collected  by  two  methods.  The  first  is  to 
place  the  source  and  receiver  on  the  seabed  and  directly  generate  and  receive  shear 
waves  as  on  land.  This  method  is  limited  to  shallow  water  (Neidell,  1986).  The  other 
method  involves  generating  compressional  waves  (P-waves)  in  the  water,  which,  at  the 
seabed,  can  convert  to  shear  waves.  The  shear  waves  may  then  be  scattered  back  to 
the  seabed  and  subsequently  converted  back  to  compressional  waves  in  the  water, 
which  are  recorded  near  the  surface.  I  will  refer  to  waves  that  are  mode-converted 
from  compressional  to  shear,  scattered,  and  then  mode-converted  back  to  compressional 


-  5  - 


as  PSSP  waves  (Figure  1.4).  P-wave  reflections  will  be  referred  to  as  PP.  One  P  for 
the  downward  leg,  and  the  second  for  the  upward. 

The  PSSP  type  of  survey  has  the  advantage  that  it  can  be  conducted  like  a 
standard  marine  seismic  survey  with  a  few  differences  mentioned  below.  This  allows 
the  data  to  be  collected  much  more  quickly  and  efficiently  than  by  the  sea  bottom 
detector  method.  One  difference  of  the  PSSP  survey  from  a  standard  marine  survey  is 
that,  because  shear  waves  travel  slower  than  P-waves,  recording  times  must  be  longer 
to  record  reflections  from  the  same  depth  as  in  a  standard  P-wave  survey.  Also,  long 
offsets  are  required  because  the  strongest  mode  conversion  occurs  for  incidence  angles 
greater  than  the  P-wave  critical  angle.  (This  will  be  shown  in  Chapter  3).  Shallow 
water  is  helpful  because  it  decreases  the  offset  required  for  reaching  the  critical  angle  of 
incidence,  but  the  water  depth  does  not  have  to  be  as  shallow  as  in  the  water  bottom 
method.  Also,  the  higher  the  P-wave  velocity  of  the  water  bottom  the  smaller  the 
critical  angle  and  thus  shorter  offsets  are  needed,  and  the  greater  the  shear  velocity  of 
the  water  bottom,  the  greater  the  energy  in  mode-converted  waves. 

The  PSSP  type  of  marine  shear  wave  survey  has  two  major  disadvantages  over 
the  sea-bottom  receiver  survey.  First,  since  the  P-  and  S-wave  energy  are  both 
recorded  as  P-waves,  the  two  must  be  separated  during  processing.  The  separation, 
based  on  move-out  velocities,  will  be  discussed  in  Chapter  3.  The  second  disadvantage 
is  that  the  polarization  of  the  S-waves  is  lost  when  they  are  converted  to  P-waves  at 
the  water  bottom.  S-wave  polarization  is  useful  for  determining  anisotropy. 

I  had  originally  intended  to  derive  an  inversion  for  PSSP  waves  and  use  it  to 
invert  the  Marathon  data.  Before  I  could  perform  an  inversion,  however,  I  needed  to 
ascertain  whether  the  data  contained  significant  mode-converted  energy.  Chapter  3 
describes  the  methods  I  used  to  search  the  data  for  PSSP  events.  I  conclude  that  the 
data  do  not  contain  significant  mode-converted  energy.  The  offsets  need  to  be  longer 
to  record  PSSP  energy.  Without  the  PSSP  waves  contained  in  the  data,  the  inversion 
for  mode-converted  energy  is  impossible. 
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2.  COMPRESSIONAL  WAVE  INVERSION 


This  chapter  will  describe  the  inversion  of  the  Marathon  data  for  compressional 
waves.  The  data  were  inverted  using  the  prestack,  common  shot,  variable  background 
inversion  method  of  Dong  (1989).  The  first  part  of  the  chapter  will  describe  the 
implementation  of  Dong’s  inversion  routine.  The  next  part  will  discuss  several 
inversions  with  different  parameters.  Only  the  kinematics,  that  is,  the  positioning  of 
the  reflectors  (Aki  and  Richards,  1980),  of  the  inversions  will  be  considered  in  this  part. 
I  will  discuss  parameter  estimation  from  the  inversion  amplitudes  (i.e.,  dynamics)  in  the 
last  part  of  the  chapter. 


2.1  Common  Shot,  c(z,z)  Inversion 


In  this  section,  I  will  briefly  describe  Dong’s  (1989)  common  shot,  c(z,z) 
implementation  of  the  general  2.5-D  inversion  formula,  equation  (1.3).  To 
accommodate  the  generality  that  c(x,z )  implies,  the  traveltime  function,  T,  and  the 
amplitude  kernel,  6,  are  computed  by  tracing  rays  from  the  output  point  to  the  source 
and  receiver  locations.  The  ray  tracing  scheme  (Docherty,  1987)  requires  constant- 
velocity  layers  bounded  by  interfaces  that  are  continuous  and  have  continues  first  and 
second  derivatives.  Also,  each  interface  must  continue  across  the  whole  line  and  its 
depth  must  be  a  single-valued  function  of  distance.  To  model  a  pinch  out,  a  layer  must 
become  very  thin,  but  maintain  finite  thickness. 


As  mentioned  in  Chapter  1,  the  amplitude  kernel,  b,  depends  on  the  source- 
receiver  geometry.  For  common-shot  data, 


y/l/o,  +  l/ar 
2(2? r)3/2  A  {x,xs)A(x,£) 


(2.1) 


(Dong,  1989).  In  equation  (2.1),  a s  and  ar  are  running  parameters  along  the  rays  for 
which  |  dx/da\  2  =  l/c2.  a  is  the  takeoff  angle  of  the  ray  at  the  receiver  location  £, 
and  its  derivative  is  the  in-plane  spreading  factor.  A  (x,x,)  and  A(x, £)  are  the  WKBJ 
amplitudes  for  rays  from  the  output  point  to  the  source  and  from  the  output  point  to 
the  receiver,  respectively.  The  WKBJ  amplitudes  involve  transmission  coefficients  and 
spreading  factors,  and  are  never  zero. 

The  first  step  of  the  inversion  is  to  Fourier  transform  the  data  by  Fast  Fourier 
Transform  (FFT).  The  w-integral  of  the  inversion  operator, 

J du>  F(w)V7w  e-,u;rlz-£)  U„(ui,£)  ,  (2.2) 
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is  the  next  step.  It  constitutes  filtering  and  an  inverse  Fourier  transform  of  the  data. 
Equation  (2.2),  also  evaluated  by  FFT,  yields  a  table  of  the  data  indexed  by  time. 
These  steps  are  independent  of  the  output  point  so  they  need  be  done  only  once,  to 
“preprocess”  the  data  prior  to  doing  the  remainder  of  the  inversion. 

The  next  step  is  the  ^-integral, 

P(x)  =  Jd£  b(z,£)up(£)  ,  (2-3) 


where  up(t,£)  is  the  preprocessed  data.  The  receiver  locations  are  discrete  so  this 
integral  is  really  a  summation  over  receiver.  The  amplitude  and  traveltime  functions 
are  determined  by  tracing  rays  from  the  output  point  to  the  source  and  to  each 
receivers.  For  each  trace,  the  data  value  indexed  by  the  traveltime  for  that  receiver  is 
linearly  interpolated  from  the  preprocessed  data  and  then  multiplied  by  the  amplitude 
function.  The  weighted  data  values  from  all  traces  are  summed  for  the  inversion 
output.  This  is  repeated  for  every  desired  output  point.  The  inversion  algorithm  may 
be  summarized  by  the  following  pseudo-code  block: 


Preprocess  data 
For  each  output  point  { 

Set  sum  to  zero 

Trace  ray  from  output  point  to  source 
For  each  receiver  { 

Trace  ray  from  output  point  to  receiver 
Get  data  value  corresponding  to  traveltime 
Weight  data  value  by  amplitude 
Add  weighted  value  to  sum 
} 

} 


This  inversion  procedure  is  applied  to  each  shot  record  separately.  For  each  shot 
record,  a  partial  image  of  the  subsurface  is  formed.  To  obtain  a  complete  image,  the 
inversion  of  each  shot  record  must  be  stacked  on  output  location.  Proper  reflector 
images  will  stack  in  phase,  but  the  impulse  response  “smiles”  and  noise  will  not.  At 
best,  the  smiles  will  overlap  slightly  and  cancel  each  other.  At  worst,  they  will  not 
overlap  but  will  be  lower  order  than  the  reflect  or  images  that  stack  constructively. 

Stacking,  however,  destroys  the  amplitude  information  in  the  inversion.  Any 
parameter  estimation  from  the  inversion  amplitudes  must  be  done  using  the  unstacked 
inverted  data. 
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2.2  Inversion  of  the  Marathon  Data 


In  this  section,  I  describe  the  application  of  Dong’s  routine  to  the  Marathon  data, 
concentrating  on  the  kinematics  of  the  inversion.  I  inverted  the  Marathon  data  twice 
with  different  parameters  and  background  velocity  models.  Finally,  I  will  compare  the 
inversion  results  to  migrations  of  the  same  data  performed  by  two  different  methods. 


2.2.1  First  Inversion 

All  migration  and  inversion  techniques  require  an  estimate  or  guess  about  the 
medium  velocity.  The  input  velocity  function  for  the  Marathon  inversion  was  nearly 
the  true  velocity  function  shown  in  Figure  1.1  with  the  interfaces  bounding  the 
constant-velocity  layers  defined  by  a  cubic  spline  fit  to  control  points  obtained  by 
digitizing  the  cross  section.  In  Figure  2.1,  the  interfaces  shown  are  denoted  by 
bandlimited  delta  functions  (i.e.,  sine  functions)  centered  at  the  proper  depths  for  each 
interface.  (Note  that  the  bandlimited  delta  function  in  Figure  2.1  is  a  function  of  depth 
rather  than  normal  distance,  as  is  the  singular  function  of  the  reflector.)  The  trace 
spacing  in  Figure  2.1  is  80  ft.  By  letting  the  second  layer  become  thin  (down  to  30  ft, 
over  the  dome),  the  pinch  out  is  modeled. 

The  input  velocity  function  is  similar  to  the  actual  model  but  not  exactly  the 
same.  The  spline  fit  causes  small  unwanted  bumps  on  the  interfaces  where  the 
control-point  spacing  is  too  fine.  The  most  notable  of  these  is  at  the  base  of  the  fault 
cutting  the  third  interface.  Also,  Marathon  had  indicated  on  the  cross-section  provided 
with  the  data  that  there  was  some  uncertainty  as  to  the  true  structure  near  the  salt 
dome  flanks. 

The  first  version  of  the  inversion  program  required  that  the  number  of  output 
traces  be  the  same  as  the  number  of  input  traces.  Therefore,  each  shot  inversion 
yielded  48  output  traces.  To  image  steeper  events  I  chose  the  output  trace  spacing  to 
be  160  ft,  twice  the  receiver  spacing.  This  gives  a  larger  output  than  input  area  so  that 
events  do  not  migrate  out  of  the  section.  The  depth  sampling  interval  was  40  ft,  with 
the  first  sample  at  zero  depth.  There  were  301  depth  samples,  so  the  last  sample  was 
at  12000  ft.  The  first  output  trace  was  located  at  the  shot  position,  and  the  other 
traces  were  to  the  right  of  the  shot.  The  farthest  trace  was,  therefore,  offset  7520  ft  to 
the  right  of  the  shot  location. 

To  save  computing  time,  rays  were  not  traced  from  every  output  point  to  every 
receiver.  Instead,  rays  were  only  traced  from  every  point  on  every  fifth  output  trace  to 
every  fifth  receiver.  The  program  linearly  interpolates  amplitude  and  traveltime 
functions  from  the  values  obtained  by  the  ray  tracing.  This  should  not  degrade  the 
results  appreciably,  because  the  traveltime  function  varies  slowly. 

Figure  2.2  is  the  inversion  of  the  shot  record  shown  in  Figure  1.3.  (Note  that 
Figure  1.3  shows  the  shot  record  after  AGC.  The  inversion  is  applied  to  the  ungained 
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record.)  The  smeared  event  at  depths  less  than  1000  ft  is  from  the  direct  arrival. 
Muting  the  direct  arrival  would  have  eliminated  this  noise.  The  events  in  Figure  2.2  at 
3000  ft,  4500  ft,  and  6300  ft  are  partial  images  of  the  first  three  interfaces.  The  event 
at  11200  ft  corresponds  to  the  model  bottom.  The  fourth  reflector,  the  sawtooth,  is 
faint  and  located  at  a  depth  of  9000  ft  and  distance  of  6000  ft. 

An  inversion  similar  to  Figure  2.2  was  obtained  for  all  290  shot  records.  Each 
yields  a  different  partial  image  of  the  subsurface.  I  sorted  the  inversions  on  the  output 
trace  location  and  stacked  them  to  form  a  full  image  of  the  subsurface,  (depths  less 
than  1300  ft  are  muted  to  eliminate  the  direct  arrival  noise).  After  stacking,  the 
output  trace  spacing  is  80  ft.  Figure  2.3,  the  ungained  stack  of  all  the  individual  shot 
inversions,  contains  several  interesting  features.  The  reflectors  stack  constructively  and 
are  imaged  well,  whereas  the  concave  upward,  migration  impulse  response  smiles  (also 
known  as  “migration  smiles”)  stack  destructively  and  are  not  noticeable.  Also,  the 
amplitude  of  the  deeper  reflectors  is  less  than  the  amplitude  of  the  shallower  reflectors. 
This  decay  in  amplitude  is  also  noticeable  in  the  shot  inversion,  Figure  2.2,  and  will  be 
discussed  in  Section  2.3. 

Stacking  the  shot  inversions  degrades  the  amplitude  information.  Parameter 
estimation  from  the  amplitudes  must  be  done  before  stacking  (Bleistein  and  Cohen, 
1989).  Because  the  amplitudes  in  Figure  2.3  have  little  meaning,  no  harm  is  done  by 
applying  gain  to  see  the  lower  reflectors  better.  Figure  2.4  is  the  stacked  section  of 
Figure  2.3  after  applying  AGC  with  a  200-sample  window.  The  gain  has  also  enhanced 
noise  in  addition  to  the  deeper  reflectors,  but  the  noise  level  is  still  relatively  low. 

Comparing  Figures  2.4  and  1.1  show  that  the  reflectors  have  been  correctly 
positioned.  Only  slight  differences  occur  between  the  cross  section  and  the  inverted 
data.  Most  of  these  are  near  the  dome  flanks  where  the  cross  section  is  in  error. 
Overall,  the  inversion  performs  well  in  locating  the  reflectors,  with  only  a  few  problems 
and  shortcomings.  The  images  of  the  steep  dome  flanks  and  fault  planes  and  the  left 
portion  of  the  leftmost  sawtooth  are  weak.  Finally,  the  model  bottom  image  has  breaks 
at  distances  8500  ft  and  14000  ft,  and  it  is  not  perfectly  flat  as  it  ought  to  be.  The 
second  inversion  was  done  to  try  to  remedy  some  of  the  shortcomings  of  this  result. 

To  understand  why  the  steep  flanks  are  imaged  weakly,  consider  an  experiment  in 
a  constant  velocity  medium  with  a  single  reflector  and  a  single  source  and  receiver. 
The  envelope  of  all  reflectors  having  the  same  reflection  time  is  the  familiar  reflection 
ellipse  with  the  source  and  receiver  at  the  foci  (Figure  2.5).  If  the  reflector  has  zero 
dip,  the  specular  point  lies  below  the  midpoint  of  the  source  and  receiver.  As  the 
reflector  dip  increases,  the  specular  point  moves  farther  up  dip  and  laterally  away  from 
the  midpoint.  Consequently,  imaging  only  in  a  region  near  the  source  and  receiver,  as 
with  this  particular  example  inversion,  discriminates  against  steep  dips.  The  flanks 
could  be  imaged  better  by  increasing  the  output  range.  Because  the  program  originally 
limited  the  number  output  traces,  the  output  range  could  only  be  increased  by 
increasing  the  output  trace  spacing.  This,  however,  would  have  led  to  aliasing  the 
steeper  events,  the  very  event  sought  after  with  the  increased  output  range.  Another 
option  would  be  to  invert  the  data  several  times  with  different  output  zones.  This  was 
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more  difficult,  however,  than  changing  the  program. 

The  weak  image  of  the  left  sawtooth  is  not  the  fault  of  the  inversion.  Instead,  it  is 
due  to  the  recording  geometry.  The  left  portion  of  the  sawtooth  was  not  illuminated 
by  the  rays  in  the  common-shot  gathers  and,  hence,  no  reflections  from  it  were 
recorded.  Therefore,  the  left  part  of  the  leftmost  sawtooth  image  is  weak. 

The  remaining  problems  are  probably  due  to  inaccuracies  of  the  background 
velocity  model.  The  breaks  in  the  bottom  image  occur  directly  beneath  sawtooth 
faults,  where  the  output  is  particularly  sensitive  to  inaccuracies  in  the  model.  The 
wobble  in  the  model  bottom  image  is  also  the  result  of  slight  inaccuracies  elsewhere  in 
the  model. 


2.2.2  Second  Inversion 

The  oackground  velocity  function  and  several  inversion  parameters  were  changed 
for  the  second  inversion  to  improve  the  image.  Figure  2.6  shows  ray  tracing  in  the  first 
background  velocity  function  from  a  subsurface  point  to  equally  spaced  points  on  the 
surface.  The  bumps  from  the  spline  fit  cause  the  local  normal  to  deviate  from  the  true 
normal.  This  causes  the  raypaths  to  be  distorted  so  that  the  traveltime  function  is  not 
smooth.  To  smooth  the  interfaces,  some  control  points  have  been  removed.  Figure  2.7 
shows  the  interfaces  for  the  second  background  velocity  function.  The  layer  velocities 
are  unchanged,  but  the  spline-  caused  bumps  on  the  interfaces  in  Figure  2.1  have  been 
removed  in  Figure  2.7.  Figure  2.8  is  a  ray  tracing  plot  from  the  same  subsurface  point 
as  in  Figure  2.6,  but  using  the  second  background  velocity  function.  Notice  that  the 
raypaths  are  not  distorted;  as  a  result,  the  traveltime  function  should  be  much 
smoother  and  thus  the  inversion  image  should  be  improved. 

Some  parameters  were  also  changed  for  the  second  inversion.  The  program  was 
altered  so  that  the  user  could  specify  the  number  of  output  traces.  For  shot  records 
with  x,  <  12000  ft,  the  second  inversion  created  300  output  traces  spaced  at  80  ft  giving 
an  image  of  the  whole  model.  To  image  the  left  edge  of  the  model  from  shot  records  on 
the  right  half  of  the  model  (i.e.,  x,  >  12000  ft),  rays  had  to  cross  the  second  interface 
twice.  The  ray  tracing  routine  did  not  allow  this.  Therefore,  for  shot  records  on  the 
right  half  of  the  model,  150  traces  were  generated  with  the  First  trace  located  at 
x  =  12000  ft.  The  inversion  of  these  shots  gave  an  image  of  the  right  half  of  the  model. 
The  increased  inversion  output  aperture  should  allow  reflectors  of  any  dip  to  be  imaged 
as  long  as  the  reflections  are  recorded  in  the  data  and  are  not  aliased. 

The  other  parameter  change  concerned  the  ray  tracing  interpolation.  Instead  of 
tracing  rays  from  every  point  on  every  Fifth  output  trace  to  every  Fifth  receiver,  rays 
were  traced  from  every  point  on  every  second  output  trace  to  every  second  receiver. 
This  was  done  to  achieve  as  accurate  an  image  as  possible.  As  a  result,  the  second 
inversion  of  all  shot  records  took  3.5  times  as  much  CPU  time  as  the  First. 


-  11  - 


Figure  2.9  is  the  inversion  of  the  shot  record  shown  in  Figure  1.3.  The  inversion 
covers  the  whole  model  although  the  imaged  portion  of  the  reflectors  is  limited  to  a 
small  zone.  The  result  is  similar  to  that  in  Figure  2.2.  The  output  outside  the  zone 
where  the  reflectors  are  imaged  consists  of  migration-smile  noise.  The  noise  is  lower 
order  so  that  it  is  not  as  visible. 

Inversions  from  nearly  all  shot  records  were  stacked  to  form  a  complete  image  of 
the  model.  The  inversions  from  a  few  shot  records  were  lost  due  to  computer  failures. 
The  stacked  image  after  AGC  is  shown  in  Figure  2.10.  Compared  to  Figure  2.4,  there 
are  some  improvements  and  also  some  disappointments.  As  expected,  the  steep  flanks 
of  the  dome  are  better  imaged.  They  appear  steeper  than  in  the  cross-section  in  Figure 
1.1.  As  mentioned  previously,  the  cross-section  is  incorrect  near  the  dome  flanks.  The 
flanks  are  probably  as  steep  as  the  inversion  depicts  them.  The  fault  plane  of  the  third 
interface  has  also  been  imaged. 

The  main  disappointment  in  this  inversion  is  the  large  degree  of  migration  smile 
noise.  Migration  smiles  occupy  more  area  in  the  shot  inversion  than  do  the  reflector 
images.  Because  some  shot  inversions  were  missing,  their  contribution  to  the 
destructive  stacking  of  the  smiles  was  lost  and,  consequently,  the  smiles  are  not 
attenuated  as  much.  The  noise  has  nearly  obliterated  the  fifth  tooth  in  the  fourth 
(sawtooth)  interface.  Selective  windowing  of  the  shot  inversions  to  eliminate  smile 
noise  before  stacking  would  result  in  an  image  with  lower  noise  and  enhanced  steep 
events. 

The  first  inversion  images  the  deeper  events  better  because  the  inversion  output 
contains  fewer  smiles.  This  is  partly  a  result  of  the  narrower  output  aperture.  The 
second  inversion  images  the  steeper  events  better,  particularly  the  dome  flanks,  because 
of  the  wider  output  aperture.  The  inversion  parameters  can  be  tuned  to  give  the 
optimum  image  for  a  particular  imaging  goal  (and  exploration  target). 


2.2.3  Migrations  of  the  Marathon  Data 

In  this  section,  I  will  compare  the  kinematic  imaging  of  the  common  shot  inversion 
to  migrations  of  the  same  data  by  two  different  methods.  Figure  2.11  is  a  poststack, 
finite-difference,  time  migration  from  the  University  of  Houston.  Detailed  comparisons 
of  time  migrations  and  depth  migrations  are  not  fair  but  general  comparisons  may  have 
some  value.  Figure  2.12  is  a  prestack,  finite-difference,  depth  migration  from  Marathon 
Oil  Company.  Comparing  Dong’s  prestack  depth  inversion  to  another  prestack  depth 
migration  is  the  best  method  to  determine  its  relative  performance. 

Compare  Figure  2.11,  the  time  migration,  to  Figures  2.4  and  2.10.  Figure  2.11 
shows  much  pull  up  in  the  center  of  the  section  because  it  is  a  time  migration.  Figures 
2.4  and  2.10  do  not  show  any  pull  up  because  they  are  sufficiently  accurate  depth 
migrations.  The  pull  up  is  not  a  valid  point  of  comparison.  Also,  the  location  of 
reflectors  cannot  be  compared.  The  next  two  points,  however,  are  legitimate.  The 
fourth  and  fifth  teeth  of  the  sawtooth  interface  are  invisible  in  the  time  migration.  In 
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both  inversions,  the  fourth  sawtooth  is  imaged  while  the  fifth  is  imaged  in  Figure  2.4. 
Also,  the  time  migration  shows  two  crossing  events  at  the  peak  of  the  dome.  Figure  1.1 
shows  that  there  ought  to  be  only  one  event.  Both  inversions  show  only  one  event  at 
the  peak  of  the  dome.  For  these  two  points,  at  least,  the  inversions  image  the  model 
better  them  does  the  time  migration. 

Now  compare  Figure  2.12,  the  depth  migration,  to  Figures  2.4  and  2.10.  The 
migration  contains  more  noise  than  do  the  inversions.  Also,  the  migration  images  the 
steep  flanks  of  the  dome  poorer  than  even  the  first  inversion  (Figure  2.4).  Both 
inversions  image  the  fourth  sawtooth  well,  whereas  the  depth  migration  images  it 
poorly.  In  the  migration,  the  fifth  sawtooth  is  completely  absent.  The  first  inversion 
(Figure  2.4)  images  it  well.  In  the  second  inversion  (Figure  2.10),  it  is  obscured  by 
noise  but  still  present.  The  inversion  imaged  the  model  better  than  the  depth 
migration  did.  This  may  be  due  to  poor  parameter  choice  for  the  depth  migration 
rather  than  a  failing  of  that  migration  method. 


2.3  Parameter  Estimation  from  Inversion  Amplitudes 


Up  to  this  point,  only  the  kinematics  of  the  inversion  have  been  discussed.  The 
main  advantage  of  inversion  is  that  it  better  “corrects”  the  reflection  amplitudes. 
Equation  (1.6)  gives  the  relationship  between  the  inversion  amplitude  and  the  singularly 
dependent  reflection  coefficient.  The  acoustic  reflection  coefficient  is  given  by 
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where  c  and  C[  are  the  medium  velocities  above  and  below  the  reflector,  6  is  the 
incidence  angle,  and  6Y  is  the  acute  angle  between  the  normal  to  the  reflector  and  the 
propagation  direction  of  the  refracted  wave  in  the  lower  medium  (Bleistein,  1987). 
Solving  equation  (2.5)  for  c  {  gives: 
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The  peak  value  of  the  inversion  amplitude  is  the  product  of  the  reflection  coefficient, 
the  cosine  of  the  particular  specular  angle  of  incidence  at  the  output  point,  and  the 
area  under  the  effective  filter  of  the  data.  The  effective  filter  is  the  product  of  all 
filters  applied  to  the  data.  To  estimate  the  velocity  of  the  medium  below  an  interface, 
one  needs  to  know  the  velocity  above  the  interface,  the  specular  angle,  and  the  area 
under  the  filter.  (This  assumes  that  the  density  does  not  vary  across  the  interface.)  It 
is  typical  to  start  with  a  known  velocity  near  the  surface  and  estimate  the  velocity 
downward  across  reflectors,  working  progressively  deeper  into  the  section.  The  upper 
medium  velocity  can  be  taken  as  known.  The  specular  angle,  however,  is  not  known 
a  priori. 

Bleistein  (1987)  describes  a  method  for  determining  the  specular  angle.  Another 
inversion  operator  can  be  constructed  by  a  minor  modification  of  the  original  inversion 
operator  kernel:  the  old  amplitude  term  is  divided  by  the  sum  of  the  gradients  of 
traveltimes  to  the  source  and  receiver.  Bleistein  calls  this  new  inversion  operator  0X. 
The  peak  value  of  this  operator  is  given  by 

0X  -  R(z,0)  J F{u>)duj  .  (2.6) 


The  ratio  of  the  peak  values  of  0  and  0X  gives  the  cosine  of  the  specular  angle.  The 
inversion  routine  can  easily  be  modified  to  output  both  0  and  fly  simultaneously  at 
virtually  no  additional  cost. 

Because  of  time  constraints,  I  could  only  invert  the  one  shot  record  in  Figure  1.3 
with  both  0  and  0X .  The  output  trace  spacing  was  80  ft,  and  48  output  traces  were 
created.  It  was  not  necessary  to  have  a  large  inversion  output  zone  for  this  shot  record 
because  I  knew  from  the  earlier  inversions  (Figures  2.2  and  2.9)  where  the  reflectors 
would  be  imaged.  The  background  velocity  model  and  ray  tracing  interpolation 
parameters  were  the  same  as  for  the  second  inversion  (Section  2.2.2).  Figure  2.13  is  the 
0  inversion  of  the  shot  record.  The  0y  inversion  looks  the  same  but  is  scaled 
differently. 

First  consider  the  uppermost  reflector.  The  first  task  is  to  get  the  peak  values  on 
the  reflector,  assuming  that  the  effective  filter  of  the  data  is  a  zero-phase,  bandlimited, 
delta  function.  The  wavelet  of  the  data,  however,  is  not  zero  phase.  It  is  a  doublet, 
and  its  peak  value  is  not  obvious;  is  it  the  absolute  maximum  of  the  positive  lobe  of  the 
doublet,  or  of  the  minimum  lobe,  or  some  other  va'ue?  For  determining  the  specular 
angle,  only  the  relative,  not  absolute,  values  of  0  and  0X  matter.  After  performing  the 
ratio  of  the  0  and  0X,  only  the  cosine  remains.  However,  the  absolute  value  of  the 
inversion  peatk  will  be  an  issue  when  finding  the  lower  medium  velocity. 
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For  finding  the  specular  angle,  I  used  the  positive  peak  value.  Because  the  peak 
value  does  not  always  fall  on  a  sample  point,  I  oversampled  0  and  0X  four  times  to  get 
better  estimates.  The  oversampling  was  done  by  FFT,  padding  the  traces  with  zeros, 
and  inverse  FFT.  The  solid  line  in  Figure  2.14  shows  the  inverse  cosine  of  the  ratio  of 
the  peak  values  of  0  and  0X  for  the  first  reflector.  In  the  zone  where  the  reflector  is 
illuminated,  from  approximately  x  =  2400  ft  to  x  —  4280  ft,  this  ought  to  be  the  specular 
incidence  angle.  The  dashed  line  is  the  theoretical  specular  incidence  angle  for  the 
approximation  that  the  first  reflector  is  a  flat  plane  at  the  depth  of  2700  ft: 

0  =  tan-1  —  , 


where  z  —  2700  ft.  This  is  a  valid  approximation  for  this  shot  record  (see  Figure  1.1). 
The  maximum  error  between  the  data  incidence  angle  and  the  theoretical  incidence 
angle  is  3%  for  2640  ft<x<3600  ft.  The  error  increases  near  the  endpoints  of  the 
illuminated  zone,  as  expected.  This  method  for  determining  the  incidence  angle  gives 
good  results  for  output  points  within  the  illuminated  zone. 

Next,  I  will  describe  how  the  area  under  the  effective  filter  was  calculated.  I  will 
consider  the  effective  filter  to  have  two  parts,  the  data  filter  and  the  processing  filter. 
The  data  filter  is  the  net  filter  effect  of  the  source  signature,  the  subsurface  (e.g. 
attenuation),  the  receiver,  and  recording  filters.  The  processing  filter  is  a  trapezoidal 
bandpass  filter  with  corner  frequencies  2,  10,  50,  and  60  Hz. 

I  assumed  that  the  direct-arrival  wavelet  represents  the  data  filter.  This  implies 
that  the  attenuation  is  negligible.  To  find  the  area  under  the  data  filter,  I  stacked  the 
shot  records  from  x,  =  2000,  4000,  6000,  14000,  16000,  18000,  20000,  and  22000  ft  for 
common  offsets.  I  omitted  the  shots  near  the  middle  of  the  model  because  the  water- 
bottom  reflection  interfered  with  the  direct  arrival.  The  stack  was  normalized  by  the 
number  of  traces  stacked.  Next,  the  first  breaks  were  aligned  and  the  data  were 
truncated  after  40  samples  (160  ms).  Next,  the  traces  had  to  be  corrected  for 
spreading.  The  3-D  acoustic  Green’s  function  predicts  decay  of  0(l/r).  Least  squares 
fitting  a  curve  of  the  form  y  =  A  f~n  to  the  peak  value  of  the  amplitude  spectrum  of 
each  trace  showed  that  the  first  arrivals  decay  as  0(l/r0'8).  I  used  the  amplitude 
spectrum  because  the  time  samples  do  not  fall  on  the  same  part  of  the  wavelet  on 
every  trace.  The  frequency  samples  do  fall  on  the  same  part  of  the  spectrum,  though. 
The  least-squares  fit  for  only  the  far-offset  traces  shows  a  decay  of  nearly  0(l/r). 
This  indicates  a  near-field  effect  is  causing  the  best-fit  decay  for  the  ensemble  of  all 
traces  to  deviate  from  0(l/r).  The  piezoelectric  transducer  source  for  the  physical 
model  experiment  is  not  a  dilational,  point  source  so  inclusion  of  near  field  effects  is 
valid  even  though  water  is  an  acoustic  media.  In  the  near  field,  decay  is  0(l/r2) 
(Lighthill,  1978).  I  applied  the  divergence  correction  by  multiplying  the  traces  by 

- - — — ,  where  A  and  B  are  determined  by  least-squares  fitting  y  =  —  +  —  to 

l  +  B/(Ar)  r  rl 

the  peak  amplitude  spectrum  values.  The  ratio  B/A  =  —320.  For  r  »  fB/A\,  the 
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divergence  correction  factor  approaches  4irr,  the  standard,  far-field  correction. 

The  processing  bandpass  filter  was  applied  to  the  divergence  corrected  traces  and 
the  area  under  the  amplitude  spectrum  for  each  trace  was  computed  by  the  trapezoidal 
rule.  The  areas  were  averaged  to  obtain  the  area  under  the  effective  filter. 

Before  the  lower  velocity  can  be  estimated,  the  phase  shift  of  the  data  wavelet 
must  be  compensated  for.  One  method  to  do  this  is  to  apply  a  wavelet  shaping  filter  to 
the  data  to  convert  the  wavelet  to  zero  phase  before  inversion.  Another  is  to 
compensate  after  inversion.  I  used  a  post-inversion  method.  Let  P  be  the  peak  (i.e., 
the  positive  maximum  rather  than  the  negative  or  absolute  maximum)  of  the  direct 
arrival  wavelet  and  Pz  be  the  peak  of  the  zero-phase  equivalent  of  the  direct  arrival 
wavelet.  I  will  call  the  positive  peak,  or  the  equivalent  negative  peak  for  phase 
reversed  reflections,  the  pseudo-peak.  I  multiplied  the  pseudo- peak  value  of  the 
inversion,  0X,  by  the  ratio  Pzj P  to  give  an  approximation  of  what  the  amplitude  would 
have  been  if  the  data  wavelet  were  zero  phase.  Again,  this  implies  the  assumption  that 
the  filter  effects  of  the  subsurface  are  negligible.  The  ratio  P,/ P  was  1.4.  If  I  had 
more  time,  I  would  have  tried  the  wavelet  shaping  method  and  compared  the  results. 

We  now  have  all  the  necessary  information  to  estimate  the  lower  medium  velocity 
using  equation  (2.5) ,  where 


Here,  Ay  =  J  F{u)dui  and  0i\peak  is  the  phase-compensated  pseudo-peak  of  (3X. 

Estimating  the  second  layer  velocity  using  the  0X  values  (after  oversampling)  for  the 
first  interface  in  the  illuminated  zone  give  approximately  22000  ft/s.  This  is  an  error  of 
40  percent.  This  poor  estimate  can  be  attributed  to  at  least  two  reasons.  First, 
equation  (2.5)  is  for  an  acoustic  reflection  coefficient  when  there  is  no  density  contrast. 
The  real  model  has  density  contrasts,  and  the  seabed  layer  is  elastic.  Second,  the 
source  was  directional.  Styrofoam  cups  were  positioned  around  the  piezoelectric 
transducer  source  and  receivers  to  reduce  the  direct  arrival.  Consequently,  the 
calculated  the  area  under  the  data  filter  using  the  direct  arrival  is  too  low  causing  the 
estimated  velocity  to  be  too  high. 
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Equation  (2.7)  can  be  rearranged  to  yield: 
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Equation  (2.7)  is  in  the  form  of  a  line:  Y  —  a  +  bX,  where 

Y  —  (1  ~  -ft)2  cos20 

(1  +  /?)2  p2  c  2 


and 


and  both  X  and  Y  are  known  quantities.  X  and  Y  may  be  plotted  and  the  unknowns, 
cx  and  pi,  can  be  determined  from  the  slope  and  intercept  of  the  best  fit  line. 

Figure  2.15  shows  Y  versus  X  for  the  first  interface  and  using  the  previously 
calculated  area  under  the  filter  and  assuming  density  of  1.0  g/cm3  for  the  water.  Only 
samples  from  the  illuminated  zone  are  plotted.  The  first  sample  in  the  illuminated 
zone  is  the  end  of  the  short  “tail”.  As  the  incidence  angle  increases,  the  points  move 
clockwise.  Therefore,  the  points  near  the  end  of  the  long  tail  correspond  to  large 
incidence  angles.  Also  plotted  in  Figure  2.15,  is  the  “best  fit”  line.  It  appears  that  the 
“best  fit”  line  should  have  negative  slope  to  fit  more  points.  The  line,  however,  must 
have  a  positive  slope  for  the  lower  density  to  be  real,  and  is  defined  by  only  four  points. 
The  first,  that  is,  right  uppermost,  point  on  the  line  corresponds  to  the  first  point  in 
Figure  2.14  where  the  data  and  theoretical  incidence  angle  curves  match,  and  the  last 
point  lies  near  the  middle  of  the  range  of  good  incidence  angles.  At  least  two  reasons 
explain  why  the  points  that  define  the  “best  fit”  line  with  positive  slope  are  from  the 
left  part  of  the  illuminated  zone  where  the  incidence  angle  is  smaller.  First,  the 
acoustic  approximation  holds  better  for  small  incidence  angles.  At  larger  incidence 
angles,  more  mode-conversion  occurs  and  equation  (2.7)  deviates  more  from  the  elastic 
P-wave  reflection  coefficient.  Second,  the  recording  aperture  necessary  for  good 
amplitudes  is  smaller  for  shorter  offsets,  and  hence  smaller  incidence  angles.  The 
recording  aperture  will  be  discussed  again  later  in  this  section.  For  these  data,  the 
recording  array  was  not  long  enough  to  give  good  inversion  amplitudes  for  reflection 
points  far  from  the  source. 

The  density  and  velocity  of  the  lower  medium,  determined  from  the  slope  and 
intercept  of  the  line  in  Figure  2.15,  are  0.90  g/cm3  and  24000  ft/s.  The  velocity  is  too 
high  because  the  area  under  the  filter  is  too  low.  Figure  2.16  shows  Y  versus  X  for  the 
first  interface  and  with  the  area  under  the  filter  adjusted  so  that  the  estimated  value  of 
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the  second  velocity  is  correct.  The  estimated  lower  medium  density  is  0.85  g/cm3. 

I  next  used  the  adjusted  area  under  the  filter  and  the  estimated  density  for  the 
second  layer  to  estimate  the  velocity  and  density  of  the  third  layer  from  the  peak  of  (3\ 
on  the  second  interface.  Figure  2.17  is  the  X  versus  Y  plot.  Again,  the  line  was  fit  for 
a  positive  slope.  The  estimated  density  and  velocity  for  the  third  layer  are  0.71  g/cm3 
and  20000  ft/sec.  The  velocity  estimate  error  is  10  percent.  This  is  a  reasonable  size 
error  for  real  data. 

Figure  2.18  is  the  X  versus  Y  plot  for  the  third  interface.  The  true  third  layer 
velocity  was  used  for  computing  X  and  Y.  The  velocity  and  density  for  the  fourth 
layer  corresponding  to  the  line  on  the  figure  are  15800  ft/s  and  .97  g/cm3.  Finding  the 
best  fit  line  for  Figure  2.18  is  not  well  defined.  One  can  imagine  many  lines  fitting  the 
data  that  would  give  densities  and  velocities  with  wide  variation.  Although  the 
velocity  estimate  is  good,  the  density  estimate  is  not  consistent  with  the  density 
estimate  from  the  first  interface.  The  second  layer  and  the  fourth  layer  were 
constructed  of  the  same  material  so  the  densities  should  be  the  same.  I  used  the 
correct  velocities  for  estimating  parameters  so  that  errors  did  not  compound.  The 
density  estimates  were  used  to  estimate  the  parameters  of  the  next  layers;  I  did  not 
know  the  true  densities. 

Finally  the  fourth  interface  gave  estimates  of  the  fifth  layer  velocity  and  density 
as  16800  ft/s  and  .78  g/cm3.  I  used  the  true  velocity  and  0.85  g/cm3  for  the  density  of 
the  fourth  layer  when  estimating  the  fifth  layer  parameters. 

Velocity  estimates  from  the  inversion  of  one  shot  record  gave  reasonable  results. 
Because  the  model  densities  are  unknown,  the  validity  of  the  density  estimates  cannot 
be  verified.  Better  accuracy  could  have  been  obtained  by  using  the  inversions  from 
more  than  one  shot.  The  multiplicity  of  data  would  have  lessened  the  opportunity  of 
error  when  fitting  a  line  to  the  scatter  plots.  The  velocity  estimates  degraded  with 
depth  because  of  limited  recording  aperture  and  compounding  errors  in  density.  The 
integration  range  for  equation  (2.3)  should  ideally  run  from  — oo  to  -t-oo.  The  actual 
integration  range  is  determined  by  the  receiver  spread.  The  dominant  contribution  of 
(2.3)  is  from  the  specular  point,  with  the  endpoints  of  integration  the  next  most 
dominant  contribution.  If  the  endpoints  axe  not  fax  enough  from  the  specular  point, 
the  inversion  amplitude  is  degraded.  Cohen  (1989)  showed  empirically  that  the 
integration  range  must  increase  for  deeper  reflectors  to  maintain  good  amplitudes  for 
parameter  estimation. 

The  major  difficulty  in  estimating  parameters  was  in  finding  the  area  under  the 
data  filter.  The  area  calculated  from  the  direct  arrivals  was  poor  because  the  source 
was  directional.  Most  field  sources  are  also  directional.  For  this  data  set,  because  the 
model  velocities  were  known  beforehand,  the  area  under  the  filter  could  be 
“determined”  by  testing  several  values  and  choosing  the  one  that  gives  the  best  results. 
With  field  data,  the  velocities  are  not  known  and  this  would  not  work. 
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3.  PSSP  SEARCH 


This  chapter  will  describe  the  methods  I  used  to  search  the  Marathon  data  for 
possible  PSSP  wave  energy. 

I  limited  the  search  for  PSSP  waves  to  shot  records  with  the  shot  location,  x3,  less 
than  6000  ft.  This  is  the  left  side  of  the  model,  where  the  water  bottom  and  next 
interface  are  nearly  flat.  Over  the  salt  dome,  where  the  structure  is  more  complex,  the 
presence  of  PSSP  energy  would  be  difficult  to  verify  because  it  would  be  difficult  to 
predict  where  the  PSSP  arrivals  would  occur. 

Tatham  et  al.  (1983)  describe  a  method  of  separating  PSSP  energy  from  PP 
energy  in  the  r—p  domain.  They  slant  stack  a  shot  record  to  transform  it  to  the  r—p 
domain.  (I  discuss  slant  stacking  and  the  r—p  domain  in  Appendix  A.)  Because  the 
water  bottom  is  flat,  the  incidence  angle  at  the  water  bottom  is  equal  to  the  take-off 
angle.  The  water  velocity  is  known  so  they  can  label  the  p-axis  with  the  incidence 
angle.  At  angles  greater  than  the  PP  critical  angle  at  the  water  bottom,  no  P-waves 
should  penetrate  the  subsurface,  only  S-waves.  So  for  p  values  beyond  the  critical 
angle,  the  only  elliptical  events  the  slant  stack  should  show  are  the  water-bottom 
reflection  and  PSSP  events.  In  Tatham ’s  slant  stack  (Figure  3.1),  the  critical  angle  is 
labeled  ic,  and  the  curved  events  on  the  right  portion  of  the  figure  are  PSSP  reflections. 
The  events  labeled  A,  B,  C,  and  D  on  the  left  side  of  the  figure  are  PP  reflections.  On 
the  right  side,  the  PSSP  reflections  from  the  reflectors  corresponding  to  the  labeled  PP 
events  are  also  labeled  A,  B,  C,  and  D.  When  inverse  slant  stacking,  Tatham  et  al. 
used  only  the  traces  with  p  corresponding  to  incidence  angles  between  40  and  75 
degrees.  The  inverse  slant  stack  contained  primarily  PSSP  reflections.  With  this 
method,  Tatham  et  al.  could  easily  identify  PSSP  energy  in  a  shot  record  and  separate 
it  from  the  PP  energy. 

I  tried  to  apply  Tatham’s  method  to  the  shot  record  with  i,.  =2000  ft  (Figure  1.3). 
I  arbitrarily  chose  this  record  as  one  from  the  section  of  the  model  where  the  first  and 
second  interfaces  are  flattest.  The  direct  arrival  and  PP  reflections  from  the  first  three 
interfaces  are  the  four  distinct  events,  and  the  PP  reflection  from  the  bottom  of  the 
model  is  less  distinct  at  about  1.53  s. 

Figure  3.2  is  the  slant  stack  of  the  shot  record  shown  in  Figure  1.3.  The  event  in 
the  upper  right  corner  corresponds  to  the  direct  arrival.  The  events  on  the  left  edge  of 
the  plot  correspond  to  the  four  PP  reflections  visible  on  the  shot  record.  The  critical 
angle  occurs  at  p  =  (15750  ft/s) —  1  =  63.5  x  10-6  s/ft.  I  will  refer  to  this  ray  parameter 
as  pc.  To  see  more  events  on  the  slant  stack,  I  applied  AGC  (Figure  3.3).  More 
curved  events  are  apparent  both  for  p  greater  than  and  less  than  pc.  The  curved 
events  axe  possible  partial  ellipses  from  reflections  (see  Appendix  A).  Linear  events  are 
also  more  pronounced  after  gaining;  these  are  edge  effects  of  the  slant  stack. 

Because  the  PP  reflection  from  the  second  interface  has  zero-offset  time  1.75  s,  a 
PSSP  reflection  from  that  interface  must  arrive  later  than  1.75  s.  None  of  the  curved 
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events  for  p  >  pe  meet  this  criterion.  For  these  data,  the  offsets  are  not  great  enough 
to  record  reflections  beyond  the  critical  angle.  Because  the  water  bottom  is  2900  ft 
deep  and  the  maximum  offset  is  4560  ft,  the  maximum  incidence  angle  at  the  water 
bottom  is  38° .  The  critical  angle  is  48° .  Therefore,  there  can  be  no  PSSP  data 
recorded  beyond  the  critical  angle  and  I  cannot  use  Tatham’s  method  for  identifying  or 
separating  the  mode-converted  energy. 

There  are,  however,  some  curved  events  in  Figure  3.3  for  p  <  pc  that  are  not 
primary  PP  reflections.  These  may  be  PSSP  reflections.  All  these  possible  PSSP 
events  occur  at  p  greater  than  11.8  x  10~6  s/ft  whereas  only  primary  PP  events  occur 
at  p  less  than  11.8  x  10-6  s/ft.  To  try  to  isolate  the  possible  PSSP  events,  I  inverted 
the  slant  stack  only  for  p  greater  than  11.8  x  10~6  s/ft.  Because  PSSP  reflections 
cannot  arrive  before  PP  reflections,  I  also  muted  the  slant  stack  for  r  less  than  0.720  s 
before  inverting  it.  Figure  3.4  shows  the  inverse  slant  stack  of  Figure  3.3  only  for  p 
greater  than  11.8  x  10~6  s/ft  and  r  greater  than  0.720  s.  The  process  of  slant  stacking 
and  inverse  slant  stacking  for  specific  ray  parameters  is  equivalent  to  dip  filtering.  The 
data  in  Figure  3.3  have  been  high-pass  dip  filtered  to  create  the  data  of  Figure  3.4. 


The  dip-filtered  shot  record  contains  many  events  that  may  be  PSSP  reflections.  I 
analyzed  the  stacking  velocity  of  these  events.  Figure  3.5  displays  the  result  of  the 
stacking  velocity  analysis  of  Figure  3.4,  a  contour  plot  of  the  semblance  across  the 
record  after  NMO  correction  has  been  applied  for  each  of  many  trial  stacking  velocities. 
Peaks  in  the  contour  plot  indicate  strongly  hyperbolic  events. 

I  varied  the  trial  stacking  velocity  from  6000  ft/s  to  11500  ft/s.  If  the  stacking 
velocity  of  a  PSSP  reflection  were  less  than  6000  ft/s,  the  reflection  would  arrive  after 
the  maximum  recording  time.  The  maximum  possible  S-wave  interval  velocity  of  the 
subwater-bottom  material  is  (15750  ft/s)/\/?  =  11140  ft/s.  The  stacking  velocity  of  a 
PSSP  reflection  from  the  second  interface,  therefore,  should  not  exceed  11500  ft/s. 
Because  dip  is  minimal  in  this  portion  of  the  section,  I  ignore  dip  effects  on  the  stacking 
velocity  and  consider  the  rms  velocity  to  be  the  stacking  velocity. 

The  velocity  analysis  (Figure  3.5)  shows  one  large  peak  and  several  smaller  peaks 
which  may  correspond  to  PSSP  reflections.  The  large  peak:  occurs  at  8500  ft/s  and  1.2 
s.  Smaller  peaks  occur  at  7500  ft/s  and  1.35  s,  7700  ft/s  and  1.0  s,  and  6700  ft/s  and 
1.75  s.  Because  the  S-wave  velocities  of  the  model  are  unknown,  I  do  not  know  which 
of  the  peaks  are  in  reasonable  locations  for  PSSP  events. 


Stacking  velocity  (actually  rms  velocity)  is  given  by 

j  1/2 
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(Sengbush,  1983),  where  t ,•  is  the  two-way  time  through  the  ith  layer  and  vt  is  the 
interval  velocity  of  the  ith  layer.  Knowing  the  thickness  of  both  layers  and  the 
interval  velocity  of  the  water  layer  allows  the  stacking  velocity  and  zero-offset  time  of  a 
PSSP  event  to  be  computed  for  a  range  of  S-wave  interval  velocities  in  the  second 
medium.  On  the  contour  plot,  I  have  plotted  a  curve  which  traces  the  possible 
locations  of  peaks  that  could  correspond  to  a  PSSP  reflection  from  the  second  interface. 
A  PSSP  event  must  have  a  peak  on  or  near  the  curve.  Because  the  large  peak  at  8500 
ft/s  and  the  small  peak  at  6700  ft/s  lie  near  the  curve  they  may  correspond  to  a  PSSP 
reflection. 

Velocity  analysis  of  common  shot  gathers  is  not  reliable,  however.  Dip  causes  the 
apex  of  the  reflection  hyperbola  to  shift  from  zero  offset.  Velocity  analysis  should  be 
performed  on  common  midpoint  gathers  where  the  hyperbola  apex  is  at  zero  offset 
regardless  of  dip.  I  dip  filtered  all  the  shot  records  for  x3  <  6000  ft  and  sorted  the  data 
into  common  midpoint  gathers.  To  improve  the  signal-to-noise  ratio,  I  stacked  five 
adjacent  CMP  gathers  on  offset.  Figure  3.6  shows  velocity  analyses  from  the  stacked 
CMP  gathers  centered  at  xm  =  1000  ft  to  xm  =  7000  ft  in  increments  of  1000  ft.  The 
velocity  axes  (horizontal)  range  from  6000  ft/s  to  11500  ft/s  for  each  of  the  individual 
contour  plots,  and  the  time  axes  (vertical)  range  from  0.6  at  the  top  to  2.0  s  at  the 
bottom.  The  curve  plotted  on  each  contour  plot  is  the  same  curve  shown  in  Figure  3.5. 
Any  possible  PSSP  event  must  have  a  semblance  peak  on  or  near  this  line,  and  the 
peak  should  be  present  consistently  in  the  individual  plots.  There  are  no  peaks  that 
consistently  fall  on  the  curve  in  Figure  3.6.  There  is  a  consistent  peak  at  8000  ft/s  and 
1.2  s  on  most  of  the  panels.  This  peak  falls  near  the  curve  but  not  on  it.  I  examined 
the  possibility  of  this  peak  corresponding  to  a  PSSP  reflection. 

The  dip-filtered  data  for  xs  <  6000  ft  were  stacked  with  a  constant  velocity  of  8000 
ft  per  second.  Figure  3.7  shows  the  constant  velocity  stack  after  AGC  has  been 
applied.  There  is  an  event  at  about  1.2  s.  The  event  is  nearly  flat  except  near  the 
edges  of  the  plot.  A  PSSP  reflection  from  the  second  interface  should  rise  significantly 
on  the  right  for  two  reasons.  First,  the  interface  is  becoming  shallower  on  the  right 
side  so  the  traveltime  should  decrease.  Also,  the  second  layer,  through  which  the  S- 
waves  travel,  is  becoming  thinner.  The  SS  leg  of  the  raypath  accounts  for  most  of  the 
traveltime  because  the  S-wave  velocity  is  much  lower  than  the  water  velocity.  As  the 
second  layer  becomes  thinner,  the  traveltime  should  decrease.  The  event  at  1.2  s  is  not 
a  PSSP  reflection  from  the  second  interface.  It  may  be  a  mis-stacked  primary  or 
multiple  P-wave  reflection. 

The  Marathon  data  set  does  not  contain  any  discernible  PSSP  reflections  above 
the  noise  level.  The  primary  reason  is  the  limited  maximum  offset,  as  I  shall  now 
demonstrate.  Brekhovskikh  (1980)  derives  the  transmission  coefficient  for  a 
compressional  wave  in  an  acoustic  medium  converting  to  a  shear  wave  in  an  elastic 
medium: 


Tps  - 


2  (p/pi)  Zt  sin  27! 

Zy  cos2  27!  +  Zt  sin2  +  Z 


(3.2) 
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In  equations  (3.2)  and  (3.3),  c  and  ci  are  the  compressional  velocities  in  the  acoustic 
and  elastic  media,  respectively,  &!  is  the  shear  velocity  in  the  elastic  medium,  and  6 ,  0X, 
and  qq  are  the  propagation  angles  of  the  incident  compressional,  transmitted 
compressional,  and  transmitted  shear  waves,  respectively.  The  density  is  p  in  the 
acoustic  medium  and  pl  in  the  elastic  medium. 

Figure  3.8  is  a  plot  of  Tp$  versus  incidence  angle,  0,  for  c  =  11750  ft/s,  Cj  =  15750 
ft/s,  and  bi  =  0.3ci.  The  densities  in  the  acoustic  and  elastic  media  are  both  unity. 
The  notch  occurs  at  the  P-wave  critical  angle.  Both  peaks  increase  in  amplitude  when 
the  shear  wave  velocity  is  a  greater  fraction  of  the  P-wave  velocity  and  the  smaller 
peak  increases  with  respect  to  the  larger  one.  Clearly,  mode  conversion  is  stronger  for 
incidence  angles  greater  than  the  critical  angle. 

Brekhovskikh  also  has  the  equation  for  the  transmission  coefficient  for  a  S-wave  in 
an  elastic  medium  to  a  P-wave  in  an  acoustic  medium.  Figure  3.9  is  a  plot  of  the  S-  to 
P-wave  transmission  coefficient  as  a  function  of  the  emergence  angle  of  the  P-wave. 
The  material  parameters  are  the  same  as  for  Figure  3.8.  Again,  the  dominant 
conversion  occurs  for  incidence/emergence  angles  greater  than  the  critical  angle. 

To  record  significant  PSSP  energy  the  offsets  must  be  long  enough  to  record  PSSP 
waves  where  the  P-wave  is  incident  on  the  water  bottom  at  post-critical  angles  and  the 
S-wave  reflects  at  the  appropriate  refraction  angle,  qq.  For  the  Marathon  data,  assume 
that  the  water  bottom  is  a  flat  plane  at  depth  Z\  —  2900  ft  and  the  second  interface  is  a 
flat  plane  at  depth  z2  =4600  ft.  Also  assume  that  the  shear  wave  velocity  is  7900  ft/s, 
or  half  the  compressional  wave  velocity,  in  the  second  layer.  The  minimum  offset  to 
record  significant  PSSP  energy  would  be 

(Zr-zjmin  =  2  zi  tan  6C  +  2  z2  tan  qq 


=  8500ft  . 


This  is  nearly  twice  the  maximum  offset  in  the  Marathon  data.  The  minimum  offset 
required  is  one  third  the  total  length  of  the  model.  Such  long  offsets  would  make  it 
difficult  to  image  the  model  with  a  PSSP  survey. 

The  long  offsets  required  for  significant  mode  conversion  and  the  large  recording 
aperture  long  offsets  require  for  accuracy  of  the  inversion  amplitudes  severly  limits  the 
usefulness  and  practicality  of  PSSP-type  marine  surveys  for  both  seabed  mapping  and 
exploration. 
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CONCLUSIONS 


Common  shot,  c(x,z ),  prestack  inversion  positions  the  reflectors  of  the  Marathon 
physical  model  data  set  accurately.  The  inversion  parameters  may  be  chosen  to  image 
the  steep  flanks  of  the  shallow  dome  or  the  deeper  sawteeth  better.  The  wider 
inversion  output  better  images  the  steeper  reflectors,  such  as  the  dome  flanks,  but  the 
sawteeth  have  been  sacrificed  to  smile  noise.  The  kinematics  of  the  inversion  compare 
favorably  to  that  of  other  migrations  of  the  data.  Inversion  theory  provides  a  means 
for  computing  the  medium  parameters  from  the  inversion  amplitudes  on  reflectors. 
However,  estimating  the  velocity  of  progressively  deeper  layers  using  the  previously 
estimated  velocities  causes  compounding  errors.  This  was  avoided  by  using  the  true 
upper  velocity  at  each  interface,  and  under  this  constraint,  estimates  from  one  shot 
record  inversion  gave  reasonable  values  for  velocity.  Because  the  model  densities  are 
unknown,  the  validity  of  the  density  estimates  cannot  be  verified.  The  parameter 
estimation  could  be  improved  by  using  the  inversion  output  from  many  shot  records. 

The  source-receiver  offsets  of  the  Marathon  data  are  not  long  enough  to  record 
reflections  incident  on  the  water  bottom  at  angles  greater  than  the  critical  angle.  For 
this  reason,  the  method  used  by  Tatham  et  al.  for  identifying  mode  converted 
reflections  and  separating  them  from  the  PP  reflections  could  not  be  used  on  the 
Marathon  data.  Stacking  velocity  analysis  of  dip-filtered  CMP  gathers  from  the 
Marathon  data  showed  one  possible  peak  that  may  have  been  the  result  of  a  PSSP 
reflection.  Stacking  proved  that  the  event  could  not  have  been  a  PSSP  reflection 
because  it  did  not  show  significant  decrease  in  traveltime  towards  the  center  of  the 
model.  There  is  no  discernible  PSSP  energy  contained  in  the  data  set.  The  minimum 
offset  required  for  recording  PSSP  waves,  assuming  flat  plane  reflectors,  is  8500  feet. 
Such  long  offsets  would  not  allow  much  of  this  model  to  be  imaged  by  PSSP  reflections. 
The  potential  usefulness  of  PSSP  type  marine  shear  wave  surveys,  both  for  exploration 
and  seabed  mapping,  is  limited  by  the  necessarily  long  offsets.  Long  offsets  require 
much  greater  aperture  and  hence,  the  ability  to  accurately  estimate  parameters  will 
require  very  long  recording  arrays  and,  therefore,  entail  operational  difficulties. 
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APPENDIX  A 


This  appendix  describes  slant  stacking,  or  Radon  transform,  and  the  r-p  domain. 
Data  in  the  t—x  domain,  u(t,  x ),  axe  slant  stacked  by  applying  the  following  operator: 


CO 

*P{t,  p)  =  f  u(t+  px,  x)dz 


(A.l) 


(Claerbout,  1985).  For  discretely  sampled  data,  equation  (A.l)  amounts  to  summing 

along  lines  of  constant  slope,  p  =  — ,  in  the  t  —x  domain  for  different  values  of  p  and  r. 

ax 

r  is  the  time  where  the  lines,  t  +  pz,  intercept  the  time  axis.  The  slope,  p,  is  also  the 

ray  parameter  p  = - .  This  implies  that  if  the  surface  velocity  is  known,  a  slant 

v 

stack  can  show  the  events  corresponding  to  a  particular  take-off  angle. 

Straight  lines  in  the  t  —  x  domain  map  to  points  in  the  r-p  domain.  The  normal 
moveout  (NMO)  hyperbola, 

(*Vv*o)2  =x2  +  4z2  > 


maps  to  an  ellipse, 


T 

2z 


2 

vs\io 


(Claerbout,  1985).  When  p  equals  zero,  r  corresponds  to  the  zero-offset  time  of  the 
reflection.  When  r  equals  zero,  p  is  the  reciprocal  of  the  NMO  velocity.  Small  p 
corresponds  to  small  take-off  angles  and  flat  events  in  t-x.  Larger  p  corresponds  to 
greater  take-off  angles  and  steeper  events. 

An  alternative  to  transforming  directly  from  the  t-x  to  the  r-p  domain,  as 
equation  (A.l)  does,  can  be  done  by  to  first  Fourier  transforming  the  original  data 
(Phinney  et  al.,  1981;  Claerbout,  1985).  Let  the  two-dimensional  Fourier  transform  of 
u(f,  z)  be  U(k,  ui): 


OO  00 

U(lc,u)=f  dt  j  dis""1'*'1  u(f,  i) 


(A.2) 


In  the  Fourier  domain,  p  =  — . 

UJ 


for  k  in  equation  (A.2), 


This  is  a  line  passing  through  the  origin.  Substituting 
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oo  oo 


U(up,  u)  =  /  -J  dxe^lt-P*)  u(t,  x)  . 


0  —  oo 


(A.3) 


Since  t  =  r+  px, 


OO  OO 


U[up,u)  =  J  dretWTJ  dxu(r  +  px,x)  . 


(A.4) 


From  equation  (A.l), 


OO 

U{up,u)  =  J  dr e,UT ip(T,  p)  . 


(A.5) 


Inverting  the  Fourier  transform  of  (A.5), 


oo 


0(r,p)  =  -~  J  due  lwT  U{up,u) 


-  OO 

00  oo  oo 

=  — ■  J  due~wrf  dxe  -’“P 1  f  dt  e u  ( «,  x) 


—  oo  —  oo 


(A.6) 


Summarizing,  the  slant  stack  can  be  obtained  by  two-dimensional  Fourier  transform  of 
the  data,  extracting  a  line  of  constant  p,  and  inverse  temporal  Fourier  transform  of  the 
data  along  the  line.  I  used  this  method  for  all  slant  stacks  in  this  thesis.  The  FFT’s 
required  padding  to  a  power  of  two  in  both  time  and  space.  To  avoid  Fourier  wrap¬ 
around,  I  padded  with  an  additional  power  of  two  in  space. 

The  inverse  of  (A. 2)  is 

OO  OO 

u(x,  t)  =  -At  f  due~'wt  f  dk  e'kl  U{k,u)  .  (A.7) 

4*  -oo  -oo 


Substituting  k  =  up  in  (A.7)  yields, 


u(x,  t) 


OO 

e~'wtJ  dpe'“Pz  U{up,u)  , 

-  OO 


(A-8) 


where  |  u  |  is  the  Jacobian  of  the  coordinate  transformation.  Substituting  for  U(up,  u) 
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from  (A. 5), 


u  (z,  f) 


CO  OO  00 

— —  f  du\uj\  f  dp  e,CJpz  f  dr e1WT 4>{r,p) 

^2J-oo  -oo  o 


(A.9) 


This  is  the  inverse  of  equation  (A.6)  There  are  two  important  differences  between  (A.9) 
and  (A.6);  the  exponents  in  the  middle  integrals  are  opposite  sign  and  (A.9)  includes 
|  w  |  whereas  (A.6)  does  not. 

By  changing  the  order  of  integration  in  (A.9)  and  applying  the  convolution  and 
shift  theorems,  an  inverse  of  the  form  of  equation  (A.l)  is  obtained: 

00 

u(x  ,t)  =  -~p(t)*f  i/>(t  -  px,  p)  dp  .  (A. 10) 

Ll T  * 

-  00 


In  equation  (A. 10),  the  asterisk  means  convolution  and  p(t)  (known  as  the  rho  filter)  is 
the  inverse  Fourier  transform  of  I  u  \ . 


Because  of  the  similarities  of  (A.9)  and  (A.6)  I  could  use  the  same  program  for 
inverse  slant  stack  as  for  forward  slant  stack.  Next,  the  rho  Filter  must  be  applied. 
|  a;  |  can  be  factored  as  [— iu\  ■  [i  sgn(o’)].  Since  tsgn(u>)  is  the  Fourier  transform  of  the 


Hilbert  transform  kernel, 


o(t)*/(t)  =  ■£*{/(<)} 


(A. 11) 


where  is  the  Hilbert  transform  of  /(f).  I  applied  the  rho  Filter  by  Hilbert 

transforming  the  data  and  the  multiplying  it  by  —iui  in  the  frequency  domain. 

As  mentioned  earlier,  a  slant  stack  can  show  the  events  corresponding  to  a 
particular  take-off  angle.  In  other  words,  a  slant  stack  decomposes  a  wave  Field  to 
plane  waves.  This  is  only  true,  however,  for  line-source  data.  A  slant  stack  of  point 
source  data  is  kinematically  a  plane-wave  decomposition.  The  amplitudes  of  the  slant 
stack  are  not  correct  for  plane-wave  decomposition  though  (Yilmaz,  1987).  Several 
authors  consider  plane-wave  decomposition  for  point-source  data  if  the  medium 
possesses  radial  symmetry.  (Brysk  and  McCowan,  1986;  Cabrera  and  Levy,  1984; 
Treitel  et  al.,  1982).  The  amplitude  treatment  in  point-source,  plane-wave 
decomposition  is  unnecessary  unless  one  is  concerned  with  the  amplitudes  in  r—p  space. 
A  forward  and  inverse  slant  stack,  if  done  properly,  should  give  the  correct  amplitudes 
in  the  t—x  domain. 
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Marathon  Data  Model 
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Figure  1.1:  Cross  section  of  data  model. 


Marathon  Data  Recording  Geometr 
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Figure  1.2:  Schematic  of  recording  array  for  each  shot  in  data. 
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Figure  1.3:  Sample  shot  record  from  the  Marathon  data.  AGC  with  a  200- 
sample  (0.8  s)  window  has  been  applied. 


-  32  - 


PSSP  Geometry 


Figure  1.4:  Ray  paths  follow  by  PSSP  mode  converted  waves. 
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First  Background  Model 
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Figure  2.1:  The  interfaces  of  the  constant  velocity  layers  of  the  background 

velocity  function  for  the  first  inversion. 


Stack  of  All  Shot  Inversions 
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inversions  of  all  shot  records.  The  reflectors  are  completely 
3ut  the  sawtooth  reflector  is  hard  to  see. 


Stack  of  All  Inversions  (gained) 
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Figure  2.4:  Stacked  section  of  Figure  2.3  after  applying  AGC  with  a  200-sample 
(8000  feet)  window. 


Reflection  Ellipse 


Figure  2.5:  The  envelope  of  all  reflectors  having  the  same  reflection  time.  It  is 
an  ellipse  with  the  source  and  receiver  at  the  foci.  Steep  reflectors 
would  be  farther  than  flat  reflectors  from  the  source-receiver 
midpoint. 
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Second  Background  Model 


40 


Figure  2.7:  The  smoothed  interfaces  of  the  constant  velocity  layers  of  the 
background  velocity  function  for  the  second  inversion. 


Inversion  of  Shot  Record 
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impulse  response  smiles  cover  a  greater  portion  of  the  inversion  than 
the  images. 


Stack  of  All  Shot  Inversions  (gained) 


ot  second  shot  inversions  after  AGC.  Steeper  events  and  the 
interface  fault  plane  have  been  imaged  better  than  in  Figure 
it  there  is  more  smile  noise. 
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Figure  2.11:  Finite-difference,  time  migration  by  the  University  of  Houston.  The 
fourth  and  fifth  sawteeth  are  not  imaged. 
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Figure  2.12:  Prestack,  finite-difference,  depth  migration  by  Marathon  Oil 
Company.  The  fifth  sawtooth  is  absent. 
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Figure  2.15:  Plot  for  estimating  velocity  and  density  from  the  inversion 
amplitudes  of  the  First  reflector.  X  =  -  sin 2 0/e  2  and, 

Y=  (1  —  /?)2/(l  +  it)  2  x  cos20/p  V.  The  slope  of  the  best  Fit  line  is 
the  reciprocal  of  the  square  of  the  density.  The  ^-intercept  is  the 
reciprocal  of  the  squire  of  the  velocity. 


Figure  2.16:  Plot  for  estimating  velocity  and  density  from  the  inversion 
amplitudes  with  the  area  under  the  filter  selected  so  that  the  velocity 
is  15750  feet/s. 
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Estimating  Velocity  and  Density 
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Figure  2.17:  Plot  for  estimating  velocity  and  density  from  the  inversion 
amplitudes  of  the  second  reflector. 


Estimating  Velocity  and  Density 


Figure  2.18:  Plot  for  estimating  velocity  and  density  from  the  inversion 
amplitudes  of  the  third  reflector. 
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Figure  3.2:  Slant  stack  of  the  (ungained)  shot  record  shown  in  Figure  1.3. 
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Figure  3.3:  Slant  stack  in  Figure  3.2  after  AGC  with  a  200-sample  (0.8  s) 
window. 
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Figure  3.4:  Inverse  slant  stack  of  Figure  3.2  for  p  >  11.8  x  10  6  s/ft.  This  is  the 
shot  record  in  Figure  1.3  after  dip  filtering. 
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Figure  3.6:  Stacking  velocity  analysis  of  dip  filtered  CMP  gathers.  A  PSSP 
event  should  be  a  peak  on  or  near  the  curved  line. 
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Figure  3.6  (cont.):  Stacking  velocity  analysis  of  dip  filtered  CMP  gathers.  A 
PSSP  event  should  be  a  peak  on  or  near  the  curved  line. 
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Figure  3.7:  Constant  velocity  stack  with  vB  =  8000  feet/s.  AGC  has  been 
applied. 
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S-P  Transmission  Coefficient 


Figure  3.9:  S  to  P  transmission  coefficient  with  the  same  parameters 
3.8.  Most  conversion  occurs  for  emergence  angles  greatei 
critical  angle  of  48° . 
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